!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! This module defines various finite element spaces. 
!!
!! \n
!!
!! The finite element spaces are defined in terms of the element
!! restrictions of the basis functions (the so-called shape functions)
!! \f$\phi_i|_{\hat{K}}\f$, where \f$\{\phi_i\}\f$ denotes a generic
!! basis. This means in particular that there is no distinction
!! between continuous and discontinuous finite element spaces, since
!! continuity depends on how the element restrictions are combined
!! together, and not on the element restrictions themselves.
!!
!! As a naming convention, all the functions with name ending in
!! <tt>*_poly</tt> return \c t_sympol objects and are private to this
!! module,  while the finite element basis constructors are public and
!! return \c t_base objects.
!!
!! Some comments concerning linear functionals, degrees of
!! freedom and projections. Let us consider a (finite dimensional)
!! space \f$V\f$ with dual space \f$V^*\f$ and a basis
!! \f$\varepsilon^\alpha\f$ of \f$V^*\f$ (i.e. a collection of
!! \f${\rm dim}(V^*)={\rm dim}(V)\f$ linearly independent functionals
!! on \f$V\f$). This basis induces a dual basis \f$e_\alpha\f$ of
!! \f$V\f$ by the condition
!! \f{displaymath}{
!!  \left< \varepsilon^\alpha,e_\beta \right>_{V^*,V} =
!!  \delta^\alpha_\beta.
!! \f}
!! Given now \f$x\in V\f$, the contravariant representation of \f$x\f$
!! can be readily obtained as
!! \f{displaymath}{
!!  x = x^\alpha e_\alpha, \qquad x^\alpha=\left< \varepsilon^\alpha,x
!!  \right>_{V^*,V}.
!! \f}
!! The two basis \f$\varepsilon^\alpha\f$ and \f$e_\alpha\f$ are often
!! referred to as (generalized) degrees of freedom of the space
!! \f$V\f$. Let us now introduce a larger space \f$\tilde{V}\f$ such
!! that \f$V\subset\tilde{V}\f$, and let us assume that the
!! functionals \f$\varepsilon^\alpha\f$ can be prolonged on
!! \f$\tilde{V}\f$. The \f$\varepsilon^\alpha\f$ then define a
!! projector
!! \f{displaymath}{
!!  \begin{array}{rcrcl}
!!   \Pi_\varepsilon & : & \tilde{V} & \to     & V \\[1mm]
!!                   &   & \tilde{x} & \mapsto & \left<
!!                   \varepsilon^\alpha,\tilde{x} \right>_{V^*,V}
!!                   e_\alpha.
!!  \end{array}
!! \f}
!! Projectors of this form are often used for the \f$\mathbb{RT}\f$
!! and the \f$\mathbb{BDM}\f$ finite element spaces. For this reason,
!! we provide in this module two functions to generate the
!! corresponding bases \f$e_\alpha\f$.
!<----------------------------------------------------------------------
module mod_fe_spaces

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_linal, only: &
   mod_linal_initialized, &
   sort,    &
   invmat

 use mod_sympoly, only: &
   mod_sympoly_initialized, &
   t_sympol,     &
   new_sympoly,  &
   nll_sympoly,  &
   all_sympoly,  &
   red_sympoly,  &
   assignment(=),&
   operator(+),  &
   operator(-),  &
   operator(*),  &
   operator(**), &
   ev,           &
   pderj,        &
   me_int,       &
   var_change,   &
   show

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me, t_lagnodes

 use mod_base, only: &
   mod_base_initialized, &
   t_base,        &
   base

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_fe_spaces_constructor, &
   mod_fe_spaces_destructor,  &
   mod_fe_spaces_initialized, &
   rt,      & ! complete f.e. space RT
   ldgh,    & ! complete f.e. space LDG-H
   bdm,     & ! complete f.e. space BDM
   cg,      & ! continuous Galerking f.e. space
   rt_vect, & ! vector f.e. space RT (orthonormal)
   dg_vect, & ! vector f.e. space DG (orthonormal)
   dg_scal, & ! scalar f.e. space DG (orthonormal)
!   legendre,& ! one-dmensional Legendre polynomials
!   solen_vect,& ! solenoindal vector functions
!   lagr_scal, & ! Lagrangian nodal basis
!   tri_nodes, & ! equally spaced triangular nodes
   rt_vect_canonical_dofs, & ! canonical RT basis
   el_nodes   ! complete set of Lagrangian nodes (see cg)
!   adapt_dg_vect  ! "diagonalize" a DG vector basis

 private

!-----------------------------------------------------------------------

! Module variables

 ! public members
 logical, protected ::               &
   mod_fe_spaces_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_fe_spaces'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_fe_spaces_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.).or. &
       (mod_kinds_initialized.eqv..false.)   .or. &
       (mod_linal_initialized.eqv..false.)   .or. &
       (mod_sympoly_initialized.eqv..false.) .or. &
       (mod_master_el_initialized.eqv..false.).or. &
       (mod_base_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_fe_spaces_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_fe_spaces_initialized = .true.
 end subroutine mod_fe_spaces_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_fe_spaces_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_fe_spaces_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_fe_spaces_initialized = .false.
 end subroutine mod_fe_spaces_destructor

!-----------------------------------------------------------------------
 
 !> Complete \f$\mathbb{RT}_k\f$ basis. See \c rt_vect_poly and \c
 !! pk_poly for details. Notice that all the arguments specifying the
 !! quadrature formula are optional: suitable defaults are chosen in
 !! \c base, in \c mod_base.
 !<
 function rt(me,k,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k                   !< order
  integer, intent(in), optional :: deg, degs !< quadrature accuracy
  !> quadrature nodes and weights
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b
 
  !> \bug gfortran-bug: the correct code is \code
  !! b = base( me, o_s=rt_vect_poly(me%d,k) ,               &
  !!           p_s=pk_poly(me%d,k) , e_s=pk_poly(me%d-1,k), &
  !!           deg=deg, degs=degs , xig=xig , wg=wg ,       &
  !!           stab=stab , xigs=xigs , wgs=wgs )
  !! \endcode
  !! but this is not supported by \c gfortran: see the bug
  !! <tt>gfortran-trunk/alloc-comps.f90</tt>. We thus have to
  !! introduce some temporaries explicitly.
  !<
  type(t_sympol), allocatable :: bug1(:,:), bug2(:), bug3(:)
   allocate(bug1(me%d,me%d*nll_sympoly(me%d,k)+nll_sympoly(me%d-1,k)))
   allocate(bug2(nll_sympoly(me%d,k)))
   allocate(bug3(nll_sympoly(me%d-1,k)))
   bug1 = rt_vect_poly(me%d,k)
   bug2 = pk_poly(me%d,k)
   bug3 = pk_poly(me%d-1,k)
   b = base( me, o_s=bug1 ,                         &
             p_s=bug2 , e_s=bug3 ,                  &
             deg=deg, degs=degs , xig=xig , wg=wg , &
             stab=stab , xigs=xigs , wgs=wgs )
   deallocate(bug1,bug2,bug3)
  ! end of bug workaround

   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function rt
 
!-----------------------------------------------------------------------
 
 !> Complete \f$\mathbb{LDG-H}_k\f$ basis. See \c rt_vect_poly and \c
 !! pk_poly for details. Notice that all the arguments specifying the
 !! quadrature formula are optional: suitable defaults are chosen in
 !! \c base, in \c mod_base.
 !<
 function ldgh(me,k,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k                   !< order
  integer, intent(in), optional :: deg, degs !< quadrature accuracy
  !> quadrature nodes and weights
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b
 
   b = base( me, o_s=pk_vect_poly(me%d,k) ,               &
             p_s=pk_poly(me%d,k) , e_s=pk_poly(me%d-1,k), &
             deg=deg, degs=degs , xig=xig , wg=wg ,       &
             stab=stab , xigs=xigs , wgs=wgs )
   
   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function ldgh
 
!-----------------------------------------------------------------------
 
 !> Complete \f$\mathbb{BDM}_k\f$ basis. See \c rt_vect_poly and \c
 !! pk_poly for details. Notice that all the arguments specifying the
 !! quadrature formula are optional: suitable defaults are chosen in
 !! \c base, in \c mod_base.
 !<
 function bdm(me,k,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k                   !< order
  integer, intent(in), optional :: deg, degs !< quadrature accuracy
  !> quadrature nodes and weights
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b
 
   b = base( me , o_s=pk_vect_poly(me%d,k) ,                 &
             p_s=pk_poly(me%d,k-1) , e_s=pk_poly(me%d-1,k) , &
             deg=deg , degs=degs , xig=xig , wg=wg ,         &
             stab=stab , xigs=xigs , wgs=wgs )
   
   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function bdm
 
!-----------------------------------------------------------------------
 
 !> Build a \c t_base object with a vector RT basis. The basis has the
 !! same structure as describen in \c rt_vect_poly. The details of the
 !! quadrature formulas can be left unspecified, in which case
 !! suitable defaults are chosen in \c base, in \c mod_base.
 !<
 function rt_vect(me,k,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k
  integer, intent(in), optional :: deg, degs
  !> quadrature nodes and weights
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b

   b = base( me , o_s=rt_vect_poly(me%d,k) ,         &
             deg=deg , degs=degs , xig=xig , wg=wg , &
             stab=stab , xigs=xigs , wgs=wgs )

   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function rt_vect
 
!-----------------------------------------------------------------------
 
 !> Recombine a \f$\mathbb{RT}_k\f$ vector basis to obtain the basis
 !! \f$\hat{e}_\alpha\f$ corresponding to the canonical degrees of
 !! freedom.
 !!
 !! The canonical degrees of freedom are given by the functionals
 !! \f$\left\{ {\hat{\varepsilon}_s}^{\alpha_s} \right\} \cup
 !!    \left\{ {\hat{\varepsilon}_e}^{\alpha_e} \right\}\f$, with
 !! \f{equation}{
 !!    \begin{array}{rcl}
 !!    \displaystyle {\hat{\varepsilon}_s}^{\alpha_s}(\hat{q}) & = &
 !!    \displaystyle \int_{\hat{e}} \hat{q}\cdot
 !!    \hat{n}\,\hat{\eta}^{\alpha_s}\,d\hat{\sigma}, \\[3mm]
 !!    \displaystyle {\hat{\varepsilon}_e}^{\alpha_e}(\hat{q}) & = &
 !!    \displaystyle \int_{\hat{K}} \hat{q}
 !!    \cdot\hat{\phi}^{\alpha_e}\,d\hat{x},
 !!    \end{array}
 !! \f}
 !! where \f$\hat{\eta}^{\alpha_s}\f$ are the elements of the basis \c
 !! eps_s and \f$\hat{\phi}^{\alpha_e}\f$ are the elements of the
 !! basis \c eps_e.
 !!
 !! \note In a typical situation, the canonical degrees of freedom are
 !! mapped on a generic element \f$K\f$ transforming the canonical
 !! basis \f$\hat{e}_\alpha\f$ by means of the Piola transformation
 !! \f{equation}{
 !!   q(x(\hat{x})) = \frac{1}{{\rm det}(B)}B\,\hat{q}(\hat{x}).
 !! \f}
 !! This defines the basis \f$e_\alpha\f$, and it is then natural to
 !! ask what the dual basis \f$\varepsilon^\alpha\f$ looks like, since
 !! these functionals define a projection operator on
 !! \f$\mathbb{RT}_k(K)\f$, as discussed in the module introduction.
 !! The duality condition yields
 !! \f{equation}{
 !!  \varepsilon^\alpha(e_\beta) = \delta^\alpha_\beta = 
 !!  \hat{\varepsilon}^\alpha(\hat{e}_\beta)\, \qquad \forall
 !!  \alpha,\beta
 !! \f}
 !! so that, for a generic vector \f$q\f$, we have
 !! \f{equation}{
 !!  \varepsilon^\alpha(q) = 
 !!  {\rm det}(B)\,\hat{\varepsilon}^\alpha(B^{-1}q).
 !! \f}
 !! Summarizing, we obtain
 !! \f{equation}{
 !!    \begin{array}{rclcrcl}
 !!    \displaystyle \varepsilon_s^{\alpha_s}(q) & = &
 !!    \displaystyle \int_{e} q\cdot
 !!    n\,\eta^{\alpha_s}\,d\sigma & = &
 !!    \displaystyle \hat{\varepsilon}_s^{\alpha_s}(\hat{q}) & = &
 !!    \displaystyle \int_{\hat{e}} \hat{q}\cdot
 !!    \hat{n}\,\hat{\eta}^{\alpha_s}\,d\hat{\sigma} ,\\[3mm]
 !!    \displaystyle \varepsilon_e^{\alpha_e}(q) & = &
 !!    {\rm det}(B)\displaystyle \int_{K} (B^{-1}q)
 !!    \cdot(B^{-1}\phi^{\alpha_e})\,dx & = &
 !!    \displaystyle {\hat{\varepsilon}_e}^{\alpha_e}(\hat{q}) & = &
 !!    \displaystyle \int_{\hat{K}} \hat{q}
 !!    \cdot\hat{\phi}^{\alpha_e}\,d\hat{x}.
 !!    \end{array}
 !! \f}
 !!
 !! \warning There are no security checks on the input arguments \c
 !! eps_s and \c eps_e.
 !<
 function rt_vect_canonical_dofs(rt_old,eps_s,eps_e) result(rt_new)
  type(t_base), intent(in) :: rt_old
  type(t_sympol), intent(in) :: eps_s(:), eps_e(:,:)
  type(t_base) :: rt_new

  integer :: nk, is, i, ii, j, k, shift_is, pk
  real(wp) :: xv0(rt_old%me%d), bhat(rt_old%me%d,rt_old%me%d-1), det
  real(wp), dimension(rt_old%mk,rt_old%mk) :: ssk, aa
  type(t_sympol) :: no_s_sigma, dotp, e(rt_old%me%d,rt_old%mk)

   !1) build the system eps_i(e_j)
   nk = size(eps_s)
   ! 1.1) boundary functionals
   do is=1,rt_old%me%d+1
     shift_is = (is-1)*nk
     ! The idea is using symbolic integration on the side. The problem
     ! is that the normal flux no_s is known as a function of the d
     ! dimensional coordinate x, while we need to work on the d-1
     ! dimensional coordinate sigma. This can be done by introducing a
     ! suitable variable change, depending on the considered side.
     xv0 = rt_old%me%xv(:,rt_old%me%isv(1,is))
     do k=1,rt_old%me%d-1
       bhat(:,k) = rt_old%me%xv(:,rt_old%me%isv(k+1,is))-xv0
     enddo
     ! scale coefficient for the computation of the integrals
     det = rt_old%me%a(is)/rt_old%me%voldm1

     do j=1,rt_old%mk
       ! variable change: x = bhat*sigma + xv0
       no_s_sigma = var_change(rt_old%no_s(j,is),bhat,xv0)
       do i=1,nk ! test function loop
         ii = shift_is + i
         ssk(ii,j) = det*me_int(rt_old%me%d-1,no_s_sigma*eps_s(i))
       enddo
     enddo
   enddo
   ! 1.2) internal functionals
   if(rt_old%k.gt.0) then
     pk = size(eps_e,2)
     shift_is = (rt_old%me%d+1)*nk
     do j=1,rt_old%mk
       do i=1,pk
         ii = shift_is + i
         dotp = 0.0_wp
         do k=1,rt_old%me%d
           dotp = dotp + rt_old%o_s(k,j)*eps_e(k,i) ! \omega \cdot \phi
         enddo
         ssk(ii,j) = me_int(rt_old%me%d,dotp)
       enddo
     enddo
   endif

   !2) invert the matrix
   call invmat( ssk , aa )

   !3) compute the new basis
   do i=1,rt_old%mk
     e(:,i) = aa(1,i)*rt_old%o_s(:,1) ! first old vector
     do j=2,rt_old%mk ! following old vectors
       e(:,i) = red_sympoly( e(:,i) + aa(j,i)*rt_old%o_s(:,j) )
     enddo
   enddo

   !4) generate the new basis
   rt_new = base( rt_old%me , o_s=e ,                                  &
     deg=rt_old%deg , degs=rt_old%degs , xig=rt_old%xig , wg=rt_old%wg,&
                  stab=rt_old%stab , xigs=rt_old%xigs , wgs=rt_old%wgs )

   ! the field k must be set explicitly (see mod_base for details)
   rt_new%k = rt_old%k
 end function rt_vect_canonical_dofs
 
!-----------------------------------------------------------------------

! function adapt_dg_vect(dg_old,e_s,gradp_s,r_s) result(b)
! ! Combine the basis functions to obtain a basis corresponding to the
! ! chosen test functions (analogous to adapt_rt_vect, but uses BDM
! ! type degrees of freedom, so that rotational test functions are also
! ! required: see Brezzi Fortin page 114)
!  type(t_base), intent(in) :: dg_old
!  type(t_sympol), intent(in) :: e_s(:), gradp_s(:,:), r_s(:,:)
!  type(t_base) :: b
!
!  integer :: nk, isl, i, ii, j, shift_is, pk, rk
!  real(wp) :: a_mat(2,1), b_mat(2), sidef2
!  real(wp), dimension(dg_old%mk,dg_old%mk) :: ssk, aa
!  type(t_sympol) :: fun, dg_new(2,dg_old%mk)
! 
!   nk = size(e_s)
!
!   !1) build the system l_i(e_j), where l_i are the functionals
!   !corresponding to the test functions (see Brezzi-Fortin page 100
!   !and page 117)
!
!   ! 1.1) boundary conditions
!   do isl=1,3
!     shift_is = (isl-1)*nk
!     do i=1,nk
!       ii = shift_is + i
!       do j=1,dg_old%mk
!         ! define the integrand as a function of xi \in [-1,1]: this
!         ! depends on the details of the definition of the reference
!         ! element in mod_base
!         select case(isl)
!          case(1) ! x=(1-xi)/2 ; y=(1+xi)/2
!           a_mat(:,1) = (/ -0.5_wp ,  0.5_wp /)
!           b_mat      = (/  0.5_wp ,  0.5_wp /)
!           sidef2 = sqrt(2.0_wp)/2.0_wp ! |e|/2
!          case(2) ! x=0 ; y=(1-xi)/2
!           a_mat(:,1) = (/  0.0_wp , -0.5_wp /)
!           b_mat      = (/  0.0_wp ,  0.5_wp /)
!           sidef2 = 1.0_wp/2.0_wp ! |e|/2
!          case(3) ! x=(1+xi)/2 ; y=0
!           a_mat(:,1) = (/  0.5_wp ,  0.0_wp /)
!           b_mat      = (/  0.5_wp ,  0.0_wp /)
!           sidef2 = 1.0_wp/2.0_wp ! |e|/2
!         end select
!         fun = var_change(dg_old%no_s(j,isl),a_mat,b_mat) * e_s(i)
!         ssk(ii,j) = sidef2*refel_int(fun)
!       enddo
!     enddo
!   enddo
!
!   ! 1.2) internal conditions
!   if(dg_old%r.gt.0) then
!     pk = size(gradp_s,2)
!     shift_is = 3*nk
!     do i=2,pk
!       ii = shift_is + i - 1 ! skip the first constant function
!       do j=1,dg_old%mk
!         ssk(ii,j) = refel_int_0101( dg_old%o_s(1,j)*gradp_s(1,i) + &
!                                     dg_old%o_s(2,j)*gradp_s(2,i) )
!       enddo
!     enddo
!
!     if(dg_old%r.gt.1) then
!       rk = size(r_s,2)
!       shift_is = 3*nk + pk - 1
!       do i=1,rk
!         ii = shift_is + i
!         do j=1,dg_old%mk
!           ssk(ii,j) = refel_int_0101( dg_old%o_s(1,j)*r_s(1,i) + &
!                                       dg_old%o_s(2,j)*r_s(2,i) )
!         enddo
!       enddo
!     endif
!   endif
!
!   !2) invert the matrix
!   call invmat( transpose(ssk) , aa )
!
!   !3) compute the new basis
!   do i=1,dg_old%mk
!     dg_new(:,i) = aa(i,1)*dg_old%o_s(:,1)
!     do j=2,dg_old%mk
!       dg_new(:,i) = red_sympoly( dg_new(:,i)   &
!                      + aa(i,j)*dg_old%o_s(:,j) )
!     enddo
!   enddo
!
!   !4) generate the new basis
!   b = base(o_s=dg_new,deg2d=dg_old%deg2d,deg1d=dg_old%deg1d, &
!            xig=dg_old%xig,wg=dg_old%wg,xige=dg_old%xige,wge=dg_old%wge)
!
!   ! the field r must be set explicitly (see mod_base for details)
!   b%r = dg_old%r
! 
! end function adapt_dg_vect
 
!-----------------------------------------------------------------------
 
 !> Build a \c t_base object with a vector DG basis. The basis has the
 !! same structure as describen in \c pk_vect_poly. The details of the
 !! quadrature formulas can be left unspecified, in which case
 !! suitable defaults are chosen in \c base, in \c mod_base.
 !<
 function dg_vect(me,k,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k
  integer, intent(in), optional :: deg, degs
  !> quadrature nodes and weights
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b

   b = base( me , o_s=pk_vect_poly(me%d,k) ,         &
             deg=deg , degs=degs , xig=xig , wg=wg , &
             stab=stab , xigs=xigs , wgs=wgs )

   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function dg_vect

!-----------------------------------------------------------------------
 
 !> Build a \c t_base object with a scalar DG basis, i.e. a basis of
 !! \f$\mathbb{P}_k(\hat{K})\f$. The basis is hierarchical as
 !! discussed in \c pk_poly. The details of the quadrature formulas
 !! can be left unspecified, in which case suitable defaults are
 !! chosen in \c base, in \c mod_base.
 !<
 function dg_scal(me,k,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k                   !< polynomial order
  integer, intent(in), optional :: deg, degs !< quadrature accuracy
  !> quadrature nodes and weights
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b
 
  !> \bug gfortran-bug: the correct code is \code
  !! b = base( me , p_s=pk_poly(me%d,k) ,             &
  !!           deg=deg , degs=degs, xig=xig , wg=wg , &
  !!           stab=stab , xigs=xigs , wgs=wgs )
  !! \endcode
  !! but this is not supported by \c gfortran: see the bug
  !! <tt>gfortran-trunk/alloc-comps.f90</tt>. We thus have to
  !! introduce some temporaries explicitly.
  !<
  type(t_sympol), allocatable :: bugtemp(:)
   allocate(bugtemp(nll_sympoly(me%d,k)))
   bugtemp = pk_poly(me%d,k)
   b = base( me , p_s=bugtemp ,             &
             deg=deg , degs=degs, xig=xig , wg=wg , &
             stab=stab , xigs=xigs , wgs=wgs )
   deallocate(bugtemp)
  ! end of bug workaround

   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function dg_scal
 
!-----------------------------------------------------------------------
 
! function solen_vect(r,deg2d,deg1d,xig,wg,xige,wge) result(b)
! ! Construct a TYPE(t_base) object with a vector basis of the
! ! solenoidal, zero flux vectors
! ! \Phi_k(K) = \{ v \in P^k(K): \nabla\cdot v=0, \; v\cdot n=0 \}
! !
! ! for k \geq 2
! !
! ! These vectors are obtained as \arrow{k} \times \nabla B_k, where
! ! \arrow{k} is a unit vertical vector and B_k is a bubble function of
! ! degree k.
! !
! ! Notice that these vectors can only be transformed with the Piola
! ! transformation, which conserves fluxes and divergence.
!  integer, intent(in) :: r
!  integer, intent(in), optional :: deg2d, deg1d
!  real(wp), intent(in), optional :: xig(:,:), wg(:), xige(:), wge(:)
!  type(t_base) :: b
!
!  integer :: j
!  type(t_sympol), allocatable :: bubbles(:), o_s(:,:)
!  character(len=100) :: message
!  character(len=*), parameter :: &
!    this_fun_name = 'solen_vect'
!
!   if(r.lt.2) then
!     write(message,'(a,i5,a)') &
!       'Solenoidal vector basis not defined for r=',r,'.'
!     call error(this_fun_name,this_mod_name,message)
!   endif
!
!   !******* COMPILER BUG ********
!   ! Many compilers have problems with the reallocate on assignment:
!   ! we therefore allocate p_s. With a correct cmpiler, this
!   ! allocate statement should be simply deleted.
!   allocate(bubbles( ((r-1)*(r))/2 ))
!   !******* XXXXXXXX XXX ********
!   bubbles = bubble_scal_poly(r+1)
!
!   allocate(o_s(2,size(bubbles)))
!   ! Compute the vector basis from the bubbles
!   do j=1,size(bubbles)
!     o_s(1,j) = -pderj(bubbles(j),2)
!     o_s(2,j) =  pderj(bubbles(j),1)
!   enddo
!
!   b = base(o_s=o_s,deg2d=deg2d,deg1d=deg1d, &
!            xig=xig,wg=wg,xige=xige,wge=wge)
!
!   deallocate(o_s)
! 
!   ! the field r must be set explicitly (see mod_base for details)
!   b%r = r
! end function solen_vect
 
!-----------------------------------------------------------------------

 !> Build a \c t_base object with a scalar, nodal basis for the nodes
 !! \f$x_j\f$, and possibly a bubble component, i.e. a basis of
 !! \f$\mathbb{P}_k(\hat{K})\f$ such that
 !! \f$\phi_i(x_j)=\delta_{ij}\f$ joined, if required, with a
 !! collection of bubble functions \f$\left\{b_{d,k_b}\right\}\f$.
 !!
 !! There have to be exactly \f$dim\left( \mathbb{P}_k(\hat{K})
 !! \right)\f$ nodes \f$x_j\f$, which can be placed arbitrarily, as
 !! far as they uniquely identify the Lagrangian basis. The nodes
 !! shall be provided in an array of \c t_lagnodes objects, as
 !! discussed in \c mod_master_el. If the input argument \c b_k is
 !! present, the corresponding bubble functions are added to the FE
 !! space.
 !!
 !! The details of the quadrature formula can be left unspecified, in
 !! which case suitable defaults are chosen in \c base, in \c
 !! mod_base.
 !!
 !! \note The Lagrangian basis and the bubble functions are simply
 !! grouped together, otherwise they are treated independently. This
 !! means that the bubble functions do not vanish at the Lagrangian
 !! nodes \f$x_j\f$, and there are no checks that the resulting FE
 !! basis is really linearly independent. For instance, the user
 !! should pay attention to:
 !! \li dimension 1, <tt>k=2</tt>, <tt>b_k=0</tt> \f$\implies\f$
 !! wrong: the bubble \f$b_{1,0}\f$ is already present in
 !! \f$\mathbb{P}_2(\hat{K}_1)\f$;
 !! \li dimension 2, <tt>k=3</tt>, <tt>b_k=0</tt> \f$\implies\f$
 !! wrong: the bubble \f$b_{2,0}\f$ is already present in
 !! \f$\mathbb{P}_3(\hat{K}_2)\f$;
 !! \li dimension 3, <tt>k=4</tt>, <tt>b_k=0</tt> \f$\implies\f$
 !! wrong: the bubble \f$b_{3,0}\f$ is already present in
 !! \f$\mathbb{P}_4(\hat{K}_3)\f$;
 !! \li dimension 2, <tt>k=3</tt>, <tt>b_k=1</tt> \f$\implies\f$
 !! correct: the bubbles \f$b_{2,1}\f$ belong to
 !! \f$\mathbb{P}_4(\hat{K}_2)\f$.
 !<
 function cg(me,k,nodes,deg,xig,wg,b_k) result(b)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  integer, intent(in) :: k              !< polynomial order
  !> Lagrangian nodes: must be consistent with \c k and <tt>me\%d</tt>
  type(t_lagnodes), intent(in) :: nodes(:)
  integer, intent(in), optional :: deg  !< quadrature accuracy
  !> quadrature formula
  real(wp), intent(in), optional :: xig(:,:), wg(:)
  !> bubble index (specify this to add bubble functions to the finite
  !! element space).
  !<
  integer, intent(in), optional :: b_k(2)
  type(t_base) :: b

  integer :: i, j, i_shift
  real(wp), allocatable :: x(:,:), van(:,:), alpha(:,:)
  type(t_sympol), allocatable :: poly(:), bubbs(:), p_s_lag(:)

   !1) map all the nodes on the master element
   call el_nodes(me,nodes,x)

   !2) get a polynomial base
   allocate(poly(nll_sympoly(me%d,k)))
   poly = all_sympoly(me%d,k) ! monomial basis
   if(size(poly).ne.size(x,2)) return ! wrong number of nodes

   !3) construct the Vandermonde system
   allocate(van(size(poly),size(poly)),alpha(size(poly),size(poly)))
   do i=1,size(poly)
     van(i,:) = ev(poly(i),x)
   enddo
   deallocate(x)

   !4) compute the coefficients
   call invmat( van , alpha )

   !5) define the new symbolic basis
   allocate(p_s_lag(size(poly)))
   do i=1,size(poly)
     p_s_lag(i) = 0.0_wp
     do j=1,size(poly)
       p_s_lag(i) = p_s_lag(i) + alpha(i,j)*poly(j)
     enddo
   enddo
   deallocate(poly,van,alpha)

   !5) add the bubble functions
   if(present(b_k)) then
     allocate(bubbs(nll_sympoly(me%d,b_k(2))))
     bubbs = bubble_poly(me%d,b_k(2))
     ! We don't want the whole bubble basis, but only the highest
     ! order functions. Since the bubble basis is hierarchical, we
     ! take the last ones.
     if(b_k(1).eq.0) then
       i_shift = 0
     else
       i_shift = nll_sympoly(me%d,b_k(1)-1)
     endif
     ! Use poly as temporary while we reallocate p_s_lag to include
     ! the bubbles.
     allocate(poly(size(p_s_lag))); poly = p_s_lag
     deallocate(p_s_lag)
     allocate(p_s_lag( size(poly) + size(bubbs)-i_shift ))
     p_s_lag = (/ poly , bubbs(i_shift+1:size(bubbs)) /)
     deallocate(bubbs,poly)
   endif

   !6) define the new t_base object
   b = base( me , p_s=p_s_lag , deg=deg , xig=xig , wg=wg )
   deallocate(p_s_lag)
 
   ! the field k must be set explicitly (see mod_base for details)
   b%k = k
 end function cg

!-----------------------------------------------------------------------
 
 !> Map Lagrangian nodes on the master element.
 !!
 !! As discussed in \c t_lagnodes and \c lag_nodes, the Lagrangian
 !! nodes are usually defined in a hierarchical manner for all the
 !! sub-simplexes composing the master element. Sometimes, however,
 !! one needs the complete set of nodes of the master element: this
 !! function maps the nodes from the sub-simplexes to the master
 !! element.
 pure subroutine el_nodes(me,nodes,x)
  type(t_me), intent(in) :: me          !< master element \f$\hat{K}\f$
  !> Lagrangian nodes: must be consistent with <tt>me\%d</tt>
  type(t_lagnodes), intent(in) :: nodes(:)
  real(wp), allocatable, intent(out) :: x(:,:) !< Lag. nodes
 
  integer :: i, e, ifl, i_shift, ndofsd_el(0:me%d)
  integer, allocatable :: ii(:)
  real(wp), allocatable :: m(:,:)

   ! count the nodes
   ndofsd_el(0) = me%d+1          ! vertices
   do e=1,me%d-1                  ! faces
     ndofsd_el(e) = nodes(e)%nn * me%children(e)%ns
   enddo
   ndofsd_el(me%d) = nodes(me%d)%nn ! elements

   ! allocate the nodes
   allocate(x(me%d,sum(ndofsd_el)))

   ! map the nodes
   x(:,1:ndofsd_el(0)) = me%xv    ! vertices
   do e=1,me%d-1                  ! faces
     if(ndofsd_el(e).gt.0) then
       allocate( ii(nodes(e)%nn) , m(me%d,e) )
       ii = (/ (i, i=1,nodes(e)%nn) /)
       do ifl=1,me%children(e)%ns ! all the faces of dimension e
         i_shift = sum(ndofsd_el(0:e-1)) + (ifl-1)*nodes(e)%nn
         do i=1,e
           m(:,i) = me%xv(:,me%children(e)%s(i+1,ifl)) &
                   - me%xv(:,me%children(e)%s(1,ifl))
         enddo
         x(:,i_shift+ii) = matmul(m,nodes(e)%x)                    &
            + spread(me%xv(:,me%children(e)%s(1,ifl)),2,nodes(e)%nn)
       enddo
       deallocate(ii,m)
     endif
   enddo
   ! elements
   x(:,sum(ndofsd_el(0:me%d-1))+1:sum(ndofsd_el)) = nodes(me%d)%x
 
 end subroutine el_nodes
 
!-----------------------------------------------------------------------

 !> Build (the polynomials of) an \f$\mathbb{RT}_k\f$ vector basis in
 !! \f$d\f$ variables. \details The basis is built from a
 !! \f$\mathbb{P}_k\f$ basis \f$\{\phi_j\}_{j=1}^{p_{d,k}}\f$ as follows:
 !! \li \f$p_{d,k}\f$ functions of type \f$[ \phi_j , 0 , \ldots , 0]\f$
 !! \li \f$p_{d,k}\f$ functions of type \f$[ 0 , \phi_j , \ldots , 0]\f$
 !! \li \f$\ldots\f$
 !! \li \f$p_{d,k}\f$ functions of type \f$[ 0 , \ldots , 0 , \phi_j]\f$
 !! \li \f$p_{d-1,k}\f$ functions of type
 !! \f$\phi_{j(k)}[x_1,\ldots,x_d]\f$
 !!
 !! where \f$j(k)\f$ is the index of any \f$\phi_j\f$ of total degree
 !! \f$k\f$. The number of functions \f$\phi_{j(k)}\f$ can be deduced
 !! observing that for any polynomial in \f$d-1\f$ variables with
 !! total degree at most \f$k\f$ it's possible to choose the exponent
 !! of \f$x_d\f$ in order to obtain a polynomial of total degree
 !! exactly \f$k\f$. The scalar basis is computed with \c pk_poly,
 !! exploiting the fact that this functions computes a hierarchical
 !! basis.
 pure function rt_vect_poly(d,k) result(p)
  integer, intent(in) :: d, k
  type(t_sympol), allocatable :: p(:,:)

  integer :: np, npm1, i, exps(d,1)
  type(t_sympol), allocatable :: ps(:)

   np = nll_sympoly(d,k)
   npm1 = nll_sympoly(d-1,k)
   allocate( ps(np) , p(d,d*np+npm1) )
   ps = pk_poly(d,k)
 
   p = 0.0_wp
   exps = 0
   do i=1,d
     ! Pk components
     p(i,(i-1)*np+1:i*np) = ps
     ! "truly" RT components
     exps(i,1) = 1
     p(i,d*np+1:d*np+npm1) = red_sympoly(                  &
       new_sympoly( (/1.0_wp/) , exps ) * ps(np-npm1+1:np) )
     exps(i,1) = 0
   enddo

   deallocate(ps)
 
 end function rt_vect_poly
 
!-----------------------------------------------------------------------
 
 !> Build (the polynomials of) a \f$(\mathbb{P}_k)^d\f$ vector basis
 !! in \f$d\f$ variables. \details The basis is built from a
 !! \f$\mathbb{P}_k\f$ basis \f$\{\phi_j\}_{j=1}^{p_k}\f$ as follows:
 !! \li \f$p_{d,k}\f$ functions of type \f$[ \phi_j , 0 , \ldots , 0]\f$
 !! \li \f$p_{d,k}\f$ functions of type \f$[ 0 , \phi_j , \ldots , 0]\f$
 !! \li \f$\ldots\f$
 !! \li \f$p_{d,k}\f$ functions of type \f$[ 0 , \ldots , 0 , \phi_j]\f$
 !!
 !! The scalar basis is computed with \c pk_poly.
 !<
 pure function pk_vect_poly(d,k) result(p)
  integer, intent(in) :: d, k
  type(t_sympol), allocatable :: p(:,:)

  integer :: np, i
  type(t_sympol), allocatable :: ps(:) ! scalar basis
 
   np = nll_sympoly(d,k)
   allocate( ps(np) , p(d,d*np) )
   ps = pk_poly(d,k)

   p = 0.0_wp
   do i=1,d
     p(i,(i-1)*np+1:i*np) = ps
   enddo

   deallocate(ps)
 
 end function pk_vect_poly
 
!-----------------------------------------------------------------------

 !> Build an orthonormal hierarchical bubble basis. The fundamental
 !! bubble, which is identically zero on \f$\partial\hat{K}_d\f$, is
 !! \f{displaymath}{
 !!   b_{d,0} = \left(1-\sum_{i=1}^dx_i\right)\prod_{i=1}^d x_i
 !! \f}
 !! and has degree \f$d+1\f$. The higher order bubbles are computed as
 !! products of \f$b_{d,0}\f$ and polynomial functions of degree 1, 2
 !! and so on. The bubble functions of degree at most \f$k\f$ form a
 !! vector space.
 !!
 !! \note The order \c k is the order of the bubble, not the
 !! polynomial order. This means that <tt>k=0</tt> corresponds to the
 !! fundamental bubble \f$b_{d,0}\f$, <tt>k=1</tt> corresponds to all
 !! the bubbles obtained as products of \f$b_{d,0}\f$ with polynomials
 !! of degree at most 1 and so on. The resulting polynomial order of
 !! the bubble functions is thus <tt>d+1+k</tt>.
 !!
 !! \note The cardinality of a bubble basis of dimension \c d and
 !! order \c k is <tt>nll_sympoly(d,k)<tt>.
 !<
 pure function bubble_poly(d,k) result(p)
  integer, intent(in) :: d, k
  type(t_sympol), allocatable :: p(:)

  integer :: nb, i, j
  integer, allocatable :: temp(:), idx(:)

  ! ifort-bug: ifort-11.1.064 has problem with the automatic variable
  ! at higher optimizations level, so we better use an allocatable
  ! variable. This bug has not yet been reported.
  !type(t_sympol) :: p1(d+1) ! 1st order monomials
  type(t_sympol), allocatable :: p1(:) ! 1st order monomials
  allocate(p1(d+1))

   ! compute a monomial basis
   p1 = all_sympoly(d,1) ! monomial basis
   ! we have to sort the monomial basis in increasing order
   allocate( temp(d+1) , idx(d+1) )
   temp = p1%deg; idx = (/(i, i=1,d+1)/)
   call sort(temp,idx)
   p1 = p1(idx)
   deallocate( temp , idx )

   ! build the bubble array
   nb = nll_sympoly(d,k) ! number of bubbles
   allocate( p(nb) , temp(nb) , idx(nb) )
   p = all_sympoly(d,k)
   ! we have to sort this too
   temp = p%deg; idx = (/(i, i=1,nb)/)
   call sort(temp,idx)
   p = p(idx)
   deallocate( temp , idx )

   ! build the fundamental bubble
   p(1) = 1.0_wp
   do i=1,d
     p(1) = p(1)-p1(1+i) ! build the (1-x1-x2-...-xd) term
   enddo
   do i=1,d
     p(1) = p(1)*p1(1+i) ! build (1-x1-x2-...-xd)*x1*x2*...*xd
   enddo
   deallocate(p1)

   ! build the higher order bubbles
   do i=2,nb
     p(i) = p(1)*p(i)
   enddo

   ! orthonormalization
   do i=1,nb
     do j=1,i-1
       p(i) = p(i) - me_int(d,p(i)*p(j))*p(j)
     enddo
     p(i) = red_sympoly(p(i))
     p(i) = (1.0_wp/sqrt(me_int(d,p(i)**2))) * p(i)
   enddo
 
 end function bubble_poly
 
!-----------------------------------------------------------------------
 
 !> Build a hierarchical basis of \f$\mathbb{P}_k(\mathbb{R}^d)\f$
 !! which is orthonormal on the master element \f$\hat{K}_d\f$ (see
 !! mod_master_el). The dimension of the basis can be computed with \c
 !! nll_sympoly, from \c mod_sympoly. The basis is hierarchical
 !! because it's composed by a constant, followed by terms of degree
 !! 1, followed by terms of degree 2 and so on.
 !!
 !! For \f$d=1\f$ the basis is composed of Legendre polynomials.
 !<
 pure function pk_poly(d,k) result(p)
 ! Compute a DG orthonormal basis
  integer, intent(in) :: d,k !< dimension and order
  type(t_sympol), allocatable :: p(:) !< size: \c nll_sympoly(d,k)

  integer :: np, i, j
  integer, allocatable :: temp(:), idx(:)

   ! determine the dimension
   np = nll_sympoly(d,k)
   allocate( p(np) , temp(np) , idx(np) )

   ! compute a monomial basis
   p = all_sympoly(d,k) ! monomial basis
   ! we have to sort the monomial basis in increasing order
   temp = p%deg; idx = (/(i, i=1,np)/)
   call sort(temp,idx)
   p = p(idx)
   deallocate(temp,idx)

   ! orthonormalization: notice that affine transformations preserve
   ! orthogonality (but change the norm)
   do i=1,np
     do j=1,i-1
       p(i) = p(i) - me_int(d,p(i)*p(j))*p(j)
     enddo
     p(i) = red_sympoly(p(i))
     p(i) = (1.0_wp/sqrt(me_int(d,p(i)**2))) * p(i)
   enddo

 end function pk_poly
 
!-----------------------------------------------------------------------
 
 !> This function is provided for completeness, however it's
 !! restricted to dimension 1. The Legendre polynomials are computed
 !! on [-1,1].
 !<
 pure function legendre_poly(k) result(p)
 ! Compute the Legendre polynomials
  integer, intent(in) :: k
  type(t_sympol) :: p(k+1)

  integer, parameter :: exp1(1,1) = reshape((/0/),(/1,1/))
  integer, parameter :: expx(1,1) = reshape((/1/),(/1,1/))
  integer :: i
  real(wp) :: ir
  type(t_sympol) :: x
 
   p(1) = new_sympoly( (/1.0_wp/) , exp1 ) ! 1
   if(k.gt.0) then
     p(2) = new_sympoly( (/1.0_wp/) , expx ) ! x
     if(k.gt.1) then
       x = p(2) ! store the polynomial x
       do i=2,k
         ir = real(i,wp)
         p(i+1) = 1.0_wp/ir * ( (2.0_wp*ir-1.0_wp)*x*p(i) &
                               - (ir-1.0_wp)*p(i-1) )
       enddo
     endif
   endif
 
 end function legendre_poly
 
!-----------------------------------------------------------------------

 !> This subroutine tests the main functionalities of this module. It
 !! should not be used in real application, but it can be useful to
 !! check changes to the module. To use it, add it to the \c PUBLIC
 !! list and simply call it in a program like the following \code
 !! program testmod
 !!  use mod_messages
 !!  use mod_kinds
 !!  use mod_fu_manager
 !!  use mod_octave_io
 !!  use mod_linal
 !!  use mod_sympoly
 !!  use mod_octave_io_sympoly
 !!  use mod_numquad
 !!  use mod_base
 !!  use mod_fe_spaces
 !!   call mod_messages_constructor()
 !!   call mod_kinds_constructor()
 !!   call mod_fu_manager_constructor()
 !!   call mod_octave_io_constructor()
 !!   call mod_linal_constructor()
 !!   call mod_sympoly_constructor()
 !!   call mod_octave_io_sympoly_constructor()
 !!   call mod_numquad_constructor()
 !!   call mod_base_constructor()
 !!   call mod_fe_spaces_constructor()
 !!   call test_module()
 !!   call mod_fe_spaces_destructor()
 !!   call mod_base_destructor()
 !!   call mod_numquad_destructor()
 !!   call mod_octave_io_sympoly_destructor()
 !!   call mod_sympoly_destructor()
 !!   call mod_linal_destructor()
 !!   call mod_octave_io_destructor()
 !!   call mod_fu_manager_destructor()
 !!   call mod_kinds_destructor()
 !!   call mod_messages_destructor()
 !! end program testmod
 !! \endcode
 !<
 subroutine test_module()

  integer :: i, j
  type(t_sympol), allocatable :: b(:), bv(:,:)

   write(*,*) "Pk basis for d=1, k=4 (Legendre basis)"
   allocate(b(nll_sympoly(1,4)))
   b = pk_poly(1,4)
   do i=1,size(b)
     call show(b(i))
   enddo
   deallocate(b)

   write(*,*) "Pk basis for d=3, k=2"
   allocate(b(nll_sympoly(3,2)))
   b = pk_poly(3,2)
   do i=1,size(b)
     call show(b(i))
   enddo
   deallocate(b)

   write(*,*) "Bubble basis for d=1, k=1"
   allocate(b(nll_sympoly(1,3-2)))
   b = bubble_poly(1,1)
   do i=1,size(b)
     call show(b(i))
   enddo
   deallocate(b)

   write(*,*) "Bubble basis for d=3, k=1"
   allocate(b(nll_sympoly(3,5-4)))
   b = bubble_poly(3,1)
   do i=1,size(b)
     call show(b(i))
   enddo
   deallocate(b)

   write(*,*) "Pk vector basis for d=3, k=1"
   allocate(bv(3,3*nll_sympoly(3,1)))
   bv = pk_vect_poly(3,1)
   do j=1,size(bv,2)
     do i=1,size(bv,1)
       call show(bv(i,j))
     enddo
   enddo
   deallocate(bv)

   write(*,*) "RT vector basis for d=3, k=0"
   allocate(bv(3,4))
   bv = rt_vect_poly(3,0)
   do j=1,size(bv,2)
     do i=1,size(bv,1)
       call show(bv(i,j))
     enddo
   enddo
   deallocate(bv)

 end subroutine test_module

!-----------------------------------------------------------------------

end module mod_fe_spaces

